Lattice model¶
This notebook simulates the evolution of the system with local dispersal, featuring three variants of the rock-paper-scissors game. All lattice-based model were developed using the Mesa framework. With Mesa's visualization tools, we can observe the spatial structures that emerge from the interactions between the agents.
Outline
- 1. Frean and Abraham model
- 2. Model with simultaneous activation of agents
- 3. Reichenbach, Mobilia and Frey model
- 4. References
Importing the necessary libraries:
import sys
import os
import solara
import mesa
import numpy as np
import pandas as pd
from matplotlib.figure import Figure
import plotly.express as px
import plotly.figure_factory as ff
import matplotlib.pyplot as plt
from mesa.experimental import JupyterViz
sys.path.append(os.path.abspath('..'))
from lattice_models.model_rand_act import RSPRandAct
from lattice_models.model_sim_act import RSPSimAct
from lattice_models.model_mobility import RSPMobility
%matplotlib inline
/Users/irenetesta/anaconda3/envs/stochpy/lib/python3.11/site-packages/solara/validate_hooks.py:122: UserWarning: /Users/irenetesta/anaconda3/envs/stochpy/lib/python3.11/site-packages/mesa/experimental/jupyter_viz.py:163: JupyterViz: `use_state` found within a nested function created on line 155
To suppress this check, replace the line with:
grid_layout, set_grid_layout = solara.use_state(grid_layout_initial) # noqa: SH104
Make sure you understand the consequences of this, by reading about the rules of hooks at:
https://solara.dev/documentation/advanced/understanding/rules-of-hooks
warnings.warn(str(e))
1. Frean and Abraham model¶
In this model the N sites are taken to be sites in a periodic, square lattice and each agent is activated once per time-step, in random order. During each activation, an agent competes with a randomly chosen neighboring agent. If the neighboring agent can be defeated, the agent wins the competition with a probability determined by the invasion rate. If the agent wins the competition, the neighboring agent is transformed into an individual of the same species of the winning agent.
For more details, see the original paper [1].
1.1. Implementation details¶
The model is defined by the class RSPRandAct in the file ../lattice_models/model_rand_act.py:
%pycat ../lattice_models/model_rand_act.py
import mesa from .patch_rand_act import PatchRandAct from mesa.datacollection import DataCollector class RSPRandAct(mesa.Model): """ A system with three species in a competitive loop: a rock beats a pair of scissors, scissors beat a sheet of paper and paper beats a rock. """ # key: species, value: list of species that the key species can beat # e.g. ROCK (0) beats SCISSOR (1), SCISSOR (1) beats PAPER (2), PAPER (2) beats ROCK (0) rules3 = {0: [1], 1: [2], 2: [0]} rules4 = {0: [1], 1: [2], 2: [3], 3: [0]} rules5 = {0: [1,2], 1: [2,3], 2: [3,4], 3: [4,0], 4: [0,1]} def __init__(self, init0, init1, init2, init3, init4, invrate0, invrate1, invrate2, invrate3, invrate4, n_species, color_map, increase_rate=False, width=50, height=50): """ Create a new playing area of (width, height) patches. """ super().__init__() self.n_species = n_species self.color_map = color_map self.increase_rate = increase_rate if self.n_species == 3: self.init_probabilities = [init0, init1, init2] self.invasion_rates = [invrate0, invrate1, invrate2] self.rules = self.rules3 elif self.n_species == 4: self.init_probabilities = [init0, init1, init2, init3] self.invasion_rates = [invrate0, invrate1, invrate2, invrate3] self.rules = self.rules4 else: # n_species == 5 self.init_probabilities = [init0, init1, init2, init3, init4] self.invasion_rates = [invrate0, invrate1, invrate2, invrate3, invrate4] self.rules = self.rules5 # set up agent scheduler self.schedule = mesa.time.RandomActivation(self) # use a simple grid, where edges wrap around. self.grid = mesa.space.SingleGrid(width, height, torus=True) # place a patch at each location, randomly initializing it according to the given initial proportions for _, (x, y) in self.grid.coord_iter(): patch_init_state = self.random.choices(range(self.n_species), weights=self.init_probabilities, k=1)[0] patch = PatchRandAct(pos=(x, y), model=self, init_state=patch_init_state) self.grid.place_agent(patch, (x, y)) self.schedule.add(patch) self.running = True # collect data model_reporter = {} for i in range(self.n_species): model_reporter[i] = lambda model, species=i: model.count_patches(species) if self.increase_rate: model_reporter['increased_rate'] = lambda model: model.invasion_rates[0] # TODO: lambda model: model.compute_average_invrate0()= self.datacollector = DataCollector( model_reporter ) self.datacollector.collect(self) def count_patches(self, species): """ Helper method to count the number of patches in the given state. """ return sum([1 for patch in self.schedule.agents if patch.state == species]) def compute_average_invrate0(self): # TODO: delete? rates = [patch.invasion_rate for patch in self.schedule.agents if patch.state == 0] return sum(rates) / len(rates) def step(self): """ Have the scheduler advance each patch by one step. """ self.schedule.step() self.datacollector.collect(self) n_existint_species = sum([1 for i in range(self.n_species) if self.count_patches(i) == 0]) if n_existint_species == self.n_species-1: self.running = False
- The model can be simulated with 3, 4, or 5 species. When there are 4 species, species <1,3> and <2,4> do not compete with each other.
- Each agent is activated once per step, in random order, with the order reshuffled each step, through the
mesa.time.RandomActivationscheduler. - Agents are situated on a rectangular, periodic grid (i.e. torus), where each cell contains exactly at most one agent. Such a grid is created by the
mesa.space.SingleGridclass. - The grid is randomly initialized by placing agents in each cell according to the initial proportions of the species (the parameters
init*are used as weights in therandom.choicefunction).
The behavior of each agent (grid patch) is defined by the class PatchRandAct in the file ../lattice_models/patch_rand_act.py:
%pycat ../lattice_models/patch_rand_act.py
import mesa class PatchRandAct(mesa.Agent): """Represents a single patch in the simulation.""" def __init__(self, pos, model, init_state): """ Create a patch, in the given state, at the given x, y position. """ super().__init__(pos, model) self.rules = model.rules self.n_species = model.n_species self.color_map = model.color_map self.x, self.y = pos self.state = init_state self.invasion_rate = self.model.invasion_rates[self.state] # TODO: delete? self._nextState = None def step(self): """ Compute the next state of the patch according to the rules of the game and the invasion rates of the species. """ # get the neighbors neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=True) # select a random neighbor neighbor = self.random.choice(neighbors) # if the neighbor can be defeated if neighbor.state in self.rules[self.state]: # change the neighbor state according to the invasion rate win_rate = self.model.invasion_rates[self.state] new_neigh_state = self.random.choices( [self.state, neighbor.state], weights=[win_rate, 1-win_rate], k=1)[0] neighbor.state = new_neigh_state # increase species 0 invasion rate if needed if (self.model.increase_rate and self.state==0 and \ new_neigh_state==self.state and self.model.schedule.steps>100): rand = self.random.uniform(0, 1e-4) new_rate = self.model.invasion_rates[self.state] + rand if 0 < new_rate < 1: self.model.invasion_rates[self.state] = new_rate def alternative_step(self): """ Compute the next state of the patch according to the rules of the game and the invasion rates of the species. """ # get the neighbors neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=True) # select a random neighbor neighbor = self.random.choice(neighbors) # if the neighbor can be defeated if neighbor.state in self.rules[self.state]: # change the neighbor state according to the invasion rate win_rate = self.invasion_rate new_neigh_state = self.random.choices( [self.state, neighbor.state], weights=[win_rate, 1-win_rate], k=1)[0] if new_neigh_state != neighbor.state: if self.state==0: neighbor.invasion_rate = self.invasion_rate else: neighbor.invasion_rate = self.model.invasion_rates[neighbor.state] neighbor.state = new_neigh_state # increase species 0 individuals invasion rate if needed if (self.model.increase_rate and self.state==0 and \ new_neigh_state==self.state and self.model.schedule.steps>100): rand = self.random.uniform(-1e-2, 1e-2) new_rate = self.invasion_rate + rand if 0 < new_rate < 1: self.invasion_rate = new_rate
- Each patch has a Moore neighborhood, i.e. the eight cells surrounding it.
- At each step, a neighbor is randomly chosen and the agent in the patch interacts with it. According to the invasion rates, the neighboring agent may be defeated and may become the same species as the agent in the patch.
- If the parameter
increase_rateis set toTrue, after 100 steps the invasion rate of species 0 is increased whenever it replicates onto a new site. The increment is a random number uniformly distributed between 0 and 1e-4 and the new invasion rate is accepted if it is in the range $0<P_r<1$).
1.2. Simulations¶
Defining functions to run simulations and visualize the evolution of the system:
def run_model(model_class, params, steps):
model = model_class(**params)
for i in range(steps):
model.step()
return model
def agent_portrayal(cell):
'''
This function defines the appearance of the agents in the visualization.
'''
return {"color": cell.color_map[cell.state], "marker": "s", "size": 1}
def draw_grid(grid, agent_portrayal, title=None, ax=None):
'''
This function draws the lattice grid.
'''
def portray(g):
x = []
y = []
s = [] # size
c = [] # color
for i in range(g.width):
for j in range(g.height):
content = g._grid[i][j]
if not content:
continue
if not hasattr(content, "__iter__"):
# Is a single grid
content = [content]
for agent in content:
data = agent_portrayal(agent)
x.append(i)
y.append(j)
if "size" in data:
s.append(data["size"])
if "color" in data:
c.append(data["color"])
out = {"x": x, "y": y}
# This is the default value for the marker size, which auto-scales
# according to the grid area.
out["s"] = (180 / max(g.width, g.height)) ** 2
if len(s) > 0:
out["s"] = s
if len(c) > 0:
out["c"] = c
return out
if ax is None:
space_fig = Figure()
space_ax = space_fig.subplots()
else:
space_ax = ax
space_ax.set_xlim(-1, grid.width)
space_ax.set_ylim(-1, grid.height)
space_ax.scatter(**portray(grid))
space_ax.set_title(title)
space_ax.axis('off')
if ax is None:
return space_fig
def line_plot(model, title=None, ax=None):
'''
This function plots the fraction of individuals of each type over time.
'''
#all_colors = []
all_species = []
patches_types = model.n_species
if hasattr(model, 'expected_swap'):
#all_colors.append('black')
all_species.append('empty')
patches_types += 1
#all_colors.extend(['tab:red', 'purple', 'gold', 'tab:blue', 'tab:green'])
all_species.extend(['$n_r$', '$n_s$', '$n_p$', '$n_u$', '$n_v$'])
if ax is None:
fig = Figure()
axs = fig.subplots()
else:
axs = ax
model_data = model.datacollector.get_model_vars_dataframe()
tot_individuals = 0
for i in range(patches_types):
tot_individuals += model_data[i]
for i in range(patches_types):
model_data[i] = model_data[i] / tot_individuals
model_data[[i for i in range(patches_types)]].plot(
ax=axs,
style='-o',
color=[model.color_map[i] for i in range(patches_types)]
)
axs.legend([all_species[i] for i in range(patches_types)])
axs.set_xlabel('Time Step')
axs.set_ylabel('Fraction of Individuals')
axs.set_title(title)
if ax is None:
solara.FigureMatplotlib(fig)
return fig
def plot_rate_species0(model, title=None):
'''
This function plots the invasion rate of species 0 over time as well as the fraction of individuals of species 0.
'''
fig = Figure()
ax = fig.subplots()
model_data = model.datacollector.get_model_vars_dataframe()
tot_individuals = model_data[0] + model_data[1] + model_data[2]
model_data[0] = model_data[0] / tot_individuals
model_data[[0, 'increased_rate']].plot(ax=ax, style='-o', color=['tab:red', 'black'])
ax.legend(['$n_r$', 'average $P_r$'])
ax.set_xlabel('Time Step')
ax.set_title(title)
solara.FigureMatplotlib(fig)
return fig
def plot_explored_area(model, title=None):
'''
This function plots the explored area over time.
'''
fig = Figure()
ax = fig.subplots()
model_data = model.datacollector.get_model_vars_dataframe()
model_data['mean_explored_area'].plot(ax=ax, style='-o', color='black')
ax.set_xlabel('Time Step')
ax.set_ylabel('Mean Explored Area')
ax.set_title(title)
solara.FigureMatplotlib(fig)
return fig
def ternary_plot_species_evolution(df, title=None, show_markers=True):
'''
Plot the evolution of species proportions over time in a ternary plot.
Parameters
----------
df : pd.DataFrame
A DataFrame containing the species proportions over time.
title : str
The title of the plot.
'''
fig = px.scatter_ternary(
df,
a=df.columns[2],
b=df.columns[0],
c=df.columns[1],
color=df.index,
color_continuous_scale='blues',
size_max=10,
title=title
)
fig.update_layout(
width=750,
)
if show_markers:
fig.update_traces(mode='lines+markers', line=dict(color='black', width=1),
marker=dict(symbol='circle', line=dict(width=1, color='black')))
else:
fig.update_traces(mode='lines', line=dict(color='black', width=1))
fig.update_layout(coloraxis_colorbar=dict(title='Time'))
fig.show('png')
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and invasion rates $P_r=0.2, P_s=0.5, P_p=0.3$. The grid will be initialized with the sites being randomly assigned to one of the three species ($n_r=0.33, n_s=0.33, n_p=0.33$). In all the simulations, species R will be represented in red, species S in purple and species P in yellow.
color_map = {
0: 'red',
1: 'purple',
2: 'gold',
}
model_params1 = {
"height": 200,
"width": 200,
"n_species": 3,
"init0": 0.33,
"init1": 0.33,
"init2": 0.33,
"init3": 0,
"init4": 0,
"invrate0": 0.2,
"invrate1": 0.5,
"invrate2": 0.3,
"invrate3": 0,
"invrate4": 0,
"color_map": color_map
}
Running the simulation for 300 steps:
model1 = run_model(model_class=RSPRandAct, params=model_params1, steps=300)
results_df = model1.datacollector.get_model_vars_dataframe().rename(columns={0: '$n_r$', 1: '$n_s$', 2: '$n_p$'})
Defining a function to build the title for the plots using the model's parameters:
def build_title(params, steps=None):
'''
This function builds the title for the plots using the parameters of the model.
'''
title = rf"$N={params['height']} \times {params['width']},"
if 'swap_exp' in params:
title += f" empty={params['init0']:.2f}, n_r^0={params['init1']:.2f}, n_s^0={params['init2']:.2f}, n_p^0={params['init3']:.2f}"
if params['init4'] > 0:
title += f", n_u^0={params['init4']:.2f}"
if params['init5'] > 0:
title += f", n_v^0={params['init5']:.2f}"
else:
title += f" n_r^0={params['init0']:.2f}, n_s^0={params['init1']:.2f}, n_p^0={params['init2']:.2f}"
if params['init3'] > 0:
title += f", n_u^0={params['init3']:.2f}"
if params['init4'] > 0:
title += f", n_v^0={params['init4']:.2f}"
if 'invrate0' in params:
title += f", P_r={params['invrate0']}, P_s={params['invrate1']}, P_p={params['invrate2']}"
if 'invrate3' in params and params['invrate3'] > 0:
title += f", P_u={params['invrate3']}"
if 'invrate4' in params and params['invrate4'] > 0:
title += f", P_v={params['invrate4']}"
if 'swap_exp' in params:
title += f", swap\_exp={params['swap_exp']}, reproduce\_exp={params['reproduce_exp']}, select\_exp={params['select_exp']}"
if steps:
title += f", step={steps}$"
else:
title += "$"
return title
Plotting the evolution of the system:
ternary_plot_species_evolution(results_df, title=build_title(params=model_params1))
Unlike the long-range dispersal model, where the weakest species outcompetes the others in the same setting, here the dynamics stabilize, and the population densities converge towards a fixed point.
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and equal invasion rates ($Pr_=0.33, P_s=0.33, P_p=0.33$). The grid will be initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.33, n_s=0.33, n_p=0.33$).
model_params2 = {
"height": 200,
"width": 200,
"n_species": 3,
"init0": 0.33,
"init1": 0.33,
"init2": 0.33,
"init3": 0,
"init4": 0,
"invrate0": 0.33,
"invrate1": 0.33,
"invrate2": 0.33,
"invrate3": 0,
"invrate4": 0,
"color_map": color_map
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPRandAct,
model_params2,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Since the widget cannot be saved in the notebook, let's run a simulation on a square lattice with $N=200\times 200$ sites and plot the grid after 100 steps and the evolution of the species proportions over time.
model2 = run_model(model_class=RSPRandAct, params=model_params2, steps=100)
draw_grid(
model2.grid,
agent_portrayal,
title=build_title(params=model_params2, steps=100)
)
line_plot(
model2,
title=build_title(params=model_params2)
)
With equal invasion rates the populations form clumps.
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and invasion rates $P_r=0.1, P_s=0.1, P_p=0.8$. The grid is initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.1, n_s=0.8, n_p=0.1$).
model_params3 = {
"height": 200,
"width": 200,
"color_map": color_map,
"n_species": 3,
"init0": 0.1,
"init1": 0.8,
"init2": 0.1,
"init3": 0,
"init4": 0,
"invrate0": 0.1,
"invrate1": 0.1,
"invrate2": 0.8,
"invrate3": 0,
"invrate4": 0,
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPRandAct,
model_params3,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the grid after 100 steps and the evolution of the species proportions over time.
model3 = run_model(model_class=RSPRandAct, params=model_params3, steps=100)
draw_grid(
model3.grid,
agent_portrayal,
title=build_title(params=model_params3, steps=100)
)
line_plot(
model3,
title=build_title(params=model_params3)
)
When species S invasion rates is much larger than the other two, species R and P form disconnected islands amongst species S.
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites and invasion rates $P_r=0.05, P_s=0.475, P_p=0.475$. The grid is initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.475, n_s=0.475, n_p=0.05$).
model_params4 = {
"height": 200,
"width": 200,
"color_map": color_map,
"n_species": 3,
"init0": 0.475,
"init1": 0.475,
"init2": 0.05,
"init3": 0,
"init4": 0,
"invrate0": 0.05,
"invrate1": 0.475,
"invrate2": 0.475,
"invrate3": 0,
"invrate4": 0,
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPRandAct,
model_params4,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the grid after 100 steps and the evolution of the species proportions over time.
model4_100 = run_model(model_class=RSPRandAct, params=model_params4, steps=100)
draw_grid(
model4_100.grid,
agent_portrayal,
title=build_title(params=model_params4, steps=100)
)
line_plot(
model4_100,
title=build_title(params=model_params4)
)
Running a simulation and plotting the grid after 300 steps and the evolution of the species proportions over time.
model4_300 = run_model(model_class=RSPRandAct, params=model_params4, steps=300)
draw_grid(
model4_300.grid,
agent_portrayal,
title=build_title(params=model_params4, steps=300)
)
line_plot(
model4_300,
title=build_title(params=model_params4)
)
If species R invasion rate is much smaller than the other two, then species R and S have similar patchy distributions, with species P persisting in fast-moving thin fronts.
Defining the parameters to perform simulations for a range of invasion probabilities chosen to sum to unity, initializing the population densities with equal proportions:
# defining a grid of points for invasion probabilities
R, S = np.mgrid[0:1:10j, 0:1:10j]
R, S = R.ravel(), S.ravel()
P = 1 - R - S
I = np.array([R, S, P]).T
# keeping only invasion probabilities greater than 0.01
I = I[np.all(I > 0.01, axis=1)]
R = I[:, 0]
S = I[:, 1]
P = I[:, 2]
# parameters for the batch run
params_batch = {
"height": 100,
"width": 100,
"color_map": color_map,
"n_species": 3,
"init0": 0.33,
"init1": 0.33,
"init2": 0.33,
"init3": 0,
"init4": 0,
"invrate3": 0,
"invrate4": 0,
}
Setting a macro to actually run the simulations or load the results of th already executed simulations from a file (executing the simulations may take a long time):
EXECUTE = False
Performing simulations:
if EXECUTE:
# dataframe to store the results
batch_results_df = pd.DataFrame()
# iterating over the invasion probabilities
for i, (Pr, Ps, Pp) in enumerate(I):
params_batch['invrate0'] = Pr
params_batch['invrate1'] = Ps
params_batch['invrate2'] = Pp
print(f"[{i+1}/{len(I)}] Running model with parameters: {params_batch}")
# simulating with the current invasion probabilities
results = mesa.batch_run(
RSPRandAct,
params_batch,
iterations=1,
max_steps=300,
number_processes=1,
data_collection_period=1,
display_progress=True,
)
# concatenating the results of the current simulation
curr_results_df = pd.DataFrame(results)
curr_results_df['simulation_id'] = i
tot_individuals = curr_results_df[0] + curr_results_df[1] + curr_results_df[2]
curr_results_df.rename(columns={0: 'individuals_r', 1: 'individuals_s', 2: 'individuals_p'}, inplace=True)
curr_results_df['$n_r$'] = curr_results_df['individuals_r'] / tot_individuals
curr_results_df['$n_s$'] = curr_results_df['individuals_s'] / tot_individuals
curr_results_df['$n_p$'] = curr_results_df['individuals_p'] / tot_individuals
batch_results_df = pd.concat([batch_results_df, curr_results_df])
# saving the results to a csv file
batch_results_df.to_csv("../results/lattice_simulations_results.csv", index=False)
else:
batch_results_df = pd.read_csv("../results/lattice_simulations_results.csv")
Plotting a countour summarizing the results of 28 simulations that are each 300 epochs long on a square lattice with $N=100\cdot 100$ sites and showing the mean value of $n_s$ over the final 100 epochs of each model run. Data are only shown for the region that had all its invasion probabilities greater than 0.01 as outside this region, with such a small invasion probability, the grid is too small to stabilize fluctuations in the population densities.
ns_final = np.array(batch_results_df[batch_results_df['Step']>200].groupby('simulation_id')['$n_s$'].mean())
fig = ff.create_ternary_contour(
np.array([I.T[2], I.T[1], I.T[0]]),
ns_final,
pole_labels=['$P_p$', '$P_s$', '$P_r'],
interp_mode='ilr',
ncontours=15,
colorscale='Portland',
coloring=None,
showmarkers=True,
showscale=True,
title=r'$\text{Mean value of }n_s\text{ over the final 100 steps of each simulation}\\\text{(population densities initialized with equal proportions)}$',
)
fig.update_layout(width=600)
fig.show('png')
The value of $n_s$ at the mean-field fixed point varies linearly from the base of the triangle to the top. Deviations from the mean-field approximation are seen toward the corners of the triangle, but $n_s$ still rises with decreasing $P_s$.
Despite the large clumps that may develop when the invasion rates are unequal, the population densities remain close to the fixed point of the mean-field equations over much of the parameter space of the invasion probabilities and the paradox, that the most aggressive species does not have the largest population, still holds.
Let's define the parameters to perform a simulation of a system with 3 species in which the invasion rate of a species (species R) is increased over time (increase_rate=True), on a square lattice with $N=200\cdot 200$ sites and initial invasion rates $P_r=0.05, P_s=0.5, P_p=0.3$. The grid is initialized with the sites being randomly assigned to each of the three species in fixed-point proportions ($n_r=0.5, n_s=0.3, n_p=0.05$).
model_params5 = {
"height": 200,
"width": 200,
"color_map": color_map,
"n_species": 3,
"init0": 0.5,
"init1": 0.3,
"init2": 0.05,
"init3": 0,
"init4": 0,
"invrate0": 0.05,
"invrate1": 0.5,
"invrate2": 0.3,
"invrate3": 0,
"invrate4": 0,
"increase_rate": True
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPRandAct,
model_params5,
measures=[line_plot, plot_rate_species0],
agent_portrayal=agent_portrayal,
)
page
model5 = run_model(model_class=RSPRandAct, params=model_params5, steps=400)
line_plot(
model5,
title=build_title(params=model_params5)
)
plot_rate_species0(
model5,
title=build_title(params=model_params5)
)
An increase in the invasion rate of species R results in a decrease of its population.
2. Model with simultaneous activation of agents¶
This is a variant of the model where the N sites are again taken to be sites in a periodic, square lattice, but at each time step, all agents are activated simultaneously. During activation, each agent competes with all its neighboring agents. If an agent is surrounded by at least three neighbors of the species that defeat it, it is transformed into an individual of that species. Note that invasion rates are not considered in this model.
2.1. Implementation details¶
The model is defined by the class RSPSimAct in the file ../lattice_models/model_sim_act.py:
%pycat ../lattice_models/model_sim_act.py
import mesa from .patch_sim_act import PatchSimAct from mesa.datacollection import DataCollector class RSPSimAct(mesa.Model): """ A system with three species in a competitive loop: a rock beats a pair of scissors, scissors beat a sheet of paper and paper beats a rock. """ rules3 = {0: [1], 1: [2], 2: [0]} rules4 = {0: [1], 1: [2], 2: [3], 3: [0]} rules5 = {0: [1,2], 1: [2,3], 2: [3,4], 3: [4,0], 4: [0,1]} def __init__(self, init0, init1, init2, init3, init4, n_species, color_map, width=50, height=50): """ Create a new playing area of (width, height) patches. """ super().__init__() self.n_species = n_species self.color_map = color_map if self.n_species == 3: self.probabilities = [init0, init1, init2] self.threshold = 3 self.rules = self.rules3 elif self.n_species == 4: self.probabilities = [init0, init1, init2, init3] self.threshold = 2 self.rules = self.rules4 else: # n_species == 5 self.probabilities = [init0, init1, init2, init3, init4] self.threshold = 3 self.rules = self.rules5 # set up agent scheduler self.schedule = mesa.time.SimultaneousActivation(self) # use a simple grid, where edges wrap around. self.grid = mesa.space.SingleGrid(width, height, torus=True) # place a patch at each location, initializing it as ROCK, SCISSOR, or PAPER for _, (x, y) in self.grid.coord_iter(): patch_init_state = self.random.choices(range(self.n_species), weights=self.probabilities, k=1)[0] patch = PatchSimAct(pos=(x, y), model=self, init_state=patch_init_state) self.grid.place_agent(patch, (x, y)) self.schedule.add(patch) self.running = True # collect data model_reporter = {} for i in range(self.n_species): model_reporter[i] = lambda model, species=i: model.count_patches(species) self.datacollector = DataCollector( model_reporter ) self.datacollector.collect(self) def count_patches(self, species): """ Helper method to count the number of patches in the given state. """ return sum([1 for patch in self.schedule.agents if patch.state == species]) def step(self): """ Have the scheduler advance each patch by one step. """ self.schedule.step() self.datacollector.collect(self) n_existint_species = sum([1 for i in range(self.n_species) if self.count_patches(i) == 0]) if n_existint_species == self.n_species-1: self.running = False
- The model can be simulated with 3, 4, or 5 species. When there are 4 species, species <1,3> and <2,4> do not compete with each other.
- Simultaneous activation of agents is implemented through the
mesa.time.SimultaneousActivationscheduler. - Agents are situated on a rectangular, periodic grid (i.e. torus), where each cell contains exactly at most one agent. Such a grid is created by the
mesa.space.SingleGridclass. - The grid is randomly initialized by placing agents in each cell according to the initial proportions of the species (the parameters
init*are used as weights in therandom.choicefunction).
The behavior of each agent (grid patch) is defined by the class PatchSimAct in the file ../lattice_models/patch_sim_act.py:
%pycat ../lattice_models/patch_sim_act.py
import mesa import numpy as np class PatchSimAct(mesa.Agent): """Represents a single patch in the simulation.""" def __init__(self, pos, model, init_state): """ Create a patch, in the given state, at the given x, y position. """ super().__init__(pos, model) self.rules = model.rules self.n_species = model.n_species self.color_map = model.color_map self.threshold = model.threshold self.x, self.y = pos self.state = init_state self._nextState = None def step(self): """ Compute the next state of the path according to the rules of the game and the number of defeats of each type. The state is not changed here, but is just computed and stored in self._nextState, because our current state may still be necessary for our neighbors to calculate their next state. """ # assume nextState is unchanged, unless changed below self._nextState = self.state # get the neighbors neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=True) # count the number of defeats of each type defeat_counts = np.zeros(self.n_species) for neighbor in neighbors: if neighbor.state != self.state and self.state in self.rules[neighbor.state]: defeat_counts[neighbor.state] += 1 # compute the winning species winners = np.argwhere(defeat_counts == np.max(defeat_counts)).reshape(-1) self.random.shuffle(winners) for i in range(len(winners)): if defeat_counts[winners[i]] >= self.threshold: self._nextState = winners[i] break def advance(self): """ Set the state to the new state, computed in the method step(). """ self.state = self._nextState
- Each patch has a Moore neighborhood, i.e. the eight cells surrounding it.
- At each step, an agent interacts with all its neighbors. If an agent is surrounded by at least three neighbors of the species that defeat it, it is transformed into an individual of that species. In case of a tie (possible with 5 species), the new species is chosen randomly among the species that defeated the agent.
2.2. Simulations¶
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species in equal proportions (init0=0.33, init1=0.33, init2=0.33).
model_params6 = {
"height": 200,
"width": 200,
"n_species": 3,
"init0": 0.33,
"init1": 0.33,
"init2": 0.33,
"init3": 0,
"init4": 0,
"color_map": color_map
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPSimAct,
model_params6,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the grid after 200 steps and the evolution of the species proportions over time.
model6 = run_model(model_class=RSPSimAct, params=model_params6, steps=200)
draw_grid(
model6.grid,
agent_portrayal,
title=build_title(params=model_params6, steps=200)
)
line_plot(
model6,
title=build_title(params=model_params6)
)
results_df = model6.datacollector.get_model_vars_dataframe().rename(columns={0: '$n_r$', 1: '$n_s$', 2: '$n_p$'})
ternary_plot_species_evolution(results_df, title=build_title(params=model_params6))
The system dynamics reach a fixed point where the population densities are equal and it leads to the emergence of entangled rotating spiral waves.
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species according to the following probabilities init0=0.2, init1=0.4, init2=0.4.
model_params7 = {
"height": 200,
"width": 200,
"n_species": 3,
"init0": 0.2,
"init1": 0.4,
"init2": 0.4,
"init3": 0,
"init4": 0,
"color_map": color_map
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPSimAct,
model_params7,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the grid after 400 steps and the evolution of the species proportions over time.
model7 = run_model(model_class=RSPSimAct, params=model_params7, steps=400)
draw_grid(
model7.grid,
agent_portrayal,
title=build_title(params=model_params7, steps=400)
)
line_plot(
model7,
title=build_title(params=model_params7)
)
results_df = model7.datacollector.get_model_vars_dataframe().rename(columns={0: '$n_r$', 1: '$n_s$', 2: '$n_p$'})
ternary_plot_species_evolution(results_df, title=build_title(params=model_params7))
The system eventually stabilizes at a fixed point with equal population densities, resulting in the formation of intricate octagonal patterns.
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species according to the following probabilities init0=0.1, init1=0.3, init2=0.5.
model_params8 = {
"height": 200,
"width": 200,
"n_species": 3,
"init0": 0.1,
"init1": 0.3,
"init2": 0.5,
"init3": 0,
"init4": 0,
"color_map": color_map
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPSimAct,
model_params8,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the evolution of the species proportions during 100 steps.
model8 = run_model(model_class=RSPSimAct, params=model_params8, steps=100)
line_plot(
model8,
title=build_title(params=model_params8)
)
This time the species with lower initial population density (species R) dominates the system and the other two species go extinct.
Let's define the parameters to perform a simulation of a system with 4 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species in equal proportions (init0=0.25, init1=0.25, init2=0.25, init3=0.25, init4=0).
model_params9 = {
"height": 200,
"width": 200,
"n_species": 4,
"init0": 0.25,
"init1": 0.25,
"init2": 0.25,
"init3": 0.25,
"init4": 0,
"color_map": {
0: 'red',
1: 'purple',
2: 'gold',
3: 'green'
}
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPSimAct,
model_params9,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the evolution of the species proportions during 200 steps.
model9 = run_model(model_class=RSPSimAct, params=model_params9, steps=200)
draw_grid(
model9.grid,
agent_portrayal,
title=build_title(params=model_params9, steps=200)
)
line_plot(
model9,
title=build_title(params=model_params9)
)
With 4 species, the system dynamics stabilize and regular spiral patterns emerge.
Let's define the parameters to perform a simulation of a system with 5 species, on a square lattice with $N=200\cdot 200$ sites, initialized with the sites being randomly assigned to each of the three species in equal proportions (init0=0.2, init1=0.2, init2=0.2, init3=0.2, init4=0.2).
model_params10 = {
"height": 200,
"width": 200,
"n_species": 5,
"init0": 0.2,
"init1": 0.2,
"init2": 0.2,
"init3": 0.2,
"init4": 0.2,
"color_map": {
0: 'red',
1: 'purple',
2: 'gold',
3: 'green',
4: 'blue',
}
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPSimAct,
model_params10,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the evolution of the species proportions during 100 steps.
model10 = run_model(model_class=RSPSimAct, params=model_params10, steps=100)
draw_grid(
model10.grid,
agent_portrayal,
title=build_title(params=model_params10, steps=100)
)
line_plot(
model10,
title=build_title(params=model_params10)
)
With 5 species, initially distributed in equal proportions, the system oscillates and again spiral patterns emerge.
3. Reichenbach, Mobilia and Frey model¶
In this model the N sites are taken to be sites in a periodic, square lattice. Each patch can be occupied by one of the species or can be blank. Each tick, the following types of events happen at defined average rates:
- Select event: Two random neighbors compete with each other. The losing patch becomes blanks.
- Reproduce event: Two random neighbors attempt to reproduce. If one of the neighbors is blank, it acquires the color of the other. Nothing happens if neither neighbor is blank.
- Swap event: Two random neighbors swap color. This represents the organisms moving.
Therefore, this model combines the local competition and reproduction from the Frean and Abraham model with spatial migration, a ubiquitous feature of real ecosystems.
For more details, see the original paper [2].
The exact number of, for instance, swap events that occur each tick is drawn from a Poisson distribution with mean equal to $(number\_ of\_ patches) * 10 ^{swap-exponent}$. A Poisson distribution defines how many times a particular event occurs given an average rate for that event assuming that the occurrences of that event are independent. Here, the occurrences of the events are approximately independent since they're being performed by different organisms.
The events occur in a random order involving random pairs of neighbors.
See the original paper for more details: Reichenbach, Tobias, Mauro Mobilia, and Erwin Frey. "Mobility promotes and jeopardizes biodiversity in rock–paper–scissors games." Nature 448.7157 (2007): 1046-1049.
3.1. Implementation details¶
The model is defined by the class RSPMobility in the file ../lattice_models/model_mobility.py:
%pycat ../lattice_models/model_mobility.py
import mesa import numpy as np from .patch_mobility import PatchMobility from mesa.datacollection import DataCollector class RSPMobility(mesa.Model): """ Represents the 2-dimensional array of patches in Rock-Scissors-Paper Game. """ rules3 = {1: [2], 2: [3], 3: [1]} rules4 = {1: [2], 2: [3], 3: [4], 4: [1]} rules5 = {1: [2,3], 2: [3,4], 3: [4,5], 4: [5,1], 5: [1,2]} def __init__( self, init0, init1, init2, init3, init4, init5, swap_exp, reproduce_exp, select_exp, n_species, color_map, width=50, height=50 ): """ Create a new playing area of (width, height) patches. """ super().__init__() self.n_species = n_species event_repetitions = (width * height) // 3 self.expected_swap = (10**swap_exp)*event_repetitions self.expected_reproduce = (10**reproduce_exp)*event_repetitions self.expected_select = (10**select_exp)*event_repetitions self.color_map = color_map if self.n_species == 3: self.probabilities = [init0, init1, init2, init3] self.rules = self.rules3 elif self.n_species == 4: self.probabilities = [init0, init1, init2, init3, init4] self.rules = self.rules4 else: # n_species == 5 self.probabilities = [init0, init1, init2, init3, init4, init5] self.rules = self.rules5 # Set up the grid and schedule. self.schedule = mesa.time.RandomActivation(self) # Use a simple grid, where edges wrap around. self.grid = mesa.space.SingleGrid(width, height, torus=True) # Place a patch at each location, initializing it as ROCK, SCISSOR, or PAPER for _, (x, y) in self.grid.coord_iter(): patch_init_state = self.random.choices(range(0, self.n_species+1), weights=self.probabilities, k=1)[0] patch = PatchMobility(pos=(x, y), model=self, init_state=patch_init_state) self.grid.place_agent(patch, (x, y)) self.schedule.add(patch) self.running = True model_reporter = {} for i in range(self.n_species+1): model_reporter[i] = lambda model, species=i: model.count_patches(species) self.datacollector = DataCollector( model_reporter ) self.datacollector.collect(self) def count_patches(self, species): """ Helper method to count the number of patches in the given state. """ return sum([1 for patch in self.schedule.agents if patch.state == species]) def step(self): """ Compute the number of global events and have the scheduler advance each patch by one step. """ # create a list of events with random Poisson distribution self.events = [] self.events.extend(['swap'] * int(np.random.poisson(self.expected_swap))) self.events.extend(['reproduce'] * int(np.random.poisson(self.expected_reproduce))) self.events.extend(['select'] * int(np.random.poisson(self.expected_select))) # shuffle events self.random.shuffle(self.events) self.schedule.step() self.datacollector.collect(self) n_existint_species = sum([1 for i in range(1, self.n_species+1) if self.count_patches(i) == 0]) if n_existint_species == self.n_species-1: self.running = False
- The model can be simulated with 3, 4, or 5 species.
init0represents the initial proportion of empty patches. When there are 4 species, species <1,3> and <2,4> do not compete with each other. - Each agent is activated once per step, in random order, with the order reshuffled each step, through the
mesa.time.RandomActivationscheduler. - Agents are situated on a rectangular, periodic grid (i.e. torus), where each cell contains exactly at most one agent. Such a grid is created by the
mesa.space.SingleGridclass. - The grid is randomly initialized by placing agents in each cell according to the initial proportions of the species (the parameters
init*are used as weights in therandom.choicefunction). - This model uses an event-based approach. It calculates how many events of each type should occur each tick and then makes agents execute those events in a random order. Note that we have to compute the number of global events rather than the number of actions that each individual patch performs since the execution of those events has to be random between the patches. That is, a single patch can't perform all of their actions in one go: suppose they end up executing 5 swaps in one tick. The swaps would be shuffling things around locally rather than allowing for an organism to travel multiple steps. Hence, this code creates a list with an entry for each event type that should occur each tick, then shuffles that list so that the events are in a random order. We then randomly activate agents who pop an event off the list, and, together with random neighboring patches, run the corresponding event.
The behavior of each agent (grid patch) is defined by the class PatchMobility in the file ../lattice_models/patch_mobility.py:
%pycat ../lattice_models/patch_mobility.py
import mesa class PatchMobility(mesa.Agent): """Represents a single patch in the simulation.""" def __init__(self, pos, model, init_state): """ Create a patch, in the given state, at the given x, y position. """ super().__init__(pos, model) self.rules = model.rules self.n_species = model.n_species self.color_map = model.color_map self.x, self.y = pos self.state = init_state self._nextState = None def step(self): """ Compute the next state of the patch and/or of a randomly chosen neighbor, according to the event that is going to happen. """ # Get the neighbors neighbors = self.model.grid.get_neighbors((self.x, self.y), moore=False) # 4 neighbors # get the action if len(self.model.events) == 0: return action = self.model.events.pop(0) # select a random neighbor neighbor = self.random.choice(neighbors) # act if action == 'swap': self._nextState = neighbor.state neighbor.state = self.state self.state = self._nextState elif action == 'select': if neighbor.state != 0 and self.state in self.rules[neighbor.state]: self.state = 0 if self.state != 0 and neighbor.state in self.rules[self.state]: neighbor.state = 0 else: if neighbor.state == 0: neighbor.state = self.state elif self.state == 0: self.state = neighbor.state
- Each patch has a Von Neumann neighborhood, i.e. the four cells surrounding it.
- Agent state
0represents an empty patch. - At each step, a neighbor is randomly chosen and the agent in the patch interacts with it, according to the event type.
Defining a function to compute the percentage of each event type at each step:
def compute_mobility_perc(params):
total = 10**(params['swap_exp']) + 10**(params['reproduce_exp']) + 10**(params['select_exp'])
print(f"Swap %: {100*(10**(params['swap_exp'])) / total:.2f}")
print(f"Reproduce %: {100*(10**(params['reproduce_exp'])) / total:.2f}")
print(f"Select %: {100*(10**(params['select_exp'])) / total:.2f}")
3.2. Simulations¶
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the three species and empty sites in equal proportions (init0=0.25, init1=0.25, init2=0.25, init3=0.25). Black patches represent empty sites.
color_map_mobility = {
0: 'black',
1: 'red',
2: 'purple',
3: 'gold'
}
model_params_mobility = {
"height": 150,
"width": 150,
"n_species": 3,
"init0": 0.25,
"init1": 0.25,
"init2": 0.25,
"init3": 0.25,
"init4": 0,
"init5": 0,
"color_map": color_map_mobility
}
We will run simulations setting swap exponent to -1, 0, 0.5 and 1 and fixing reproduce and select exponents to 0.
Setting parameters and computing the percentage of each event type at each step.
model_params11 = model_params_mobility.copy()
model_params11['swap_exp'] = -1
model_params11['reproduce_exp'] = 0
model_params11['select_exp'] = 0
print("Configuration 1")
compute_mobility_perc(model_params11)
print("----------------")
model_params12 = model_params_mobility.copy()
model_params12['swap_exp'] = 0
model_params12['reproduce_exp'] = 0
model_params12['select_exp'] = 0
print("Configuration 2")
compute_mobility_perc(model_params12)
print("----------------")
model_params13 = model_params_mobility.copy()
model_params13['swap_exp'] = 0.5
model_params13['reproduce_exp'] = 0
model_params13['select_exp'] = 0
print("Configuration 3")
compute_mobility_perc(model_params13)
print("----------------")
model_params14 = model_params_mobility.copy()
model_params14['swap_exp'] = 1
model_params14['reproduce_exp'] = 0
model_params14['select_exp'] = 0
print("Configuration 4")
compute_mobility_perc(model_params14)
print("----------------")
Configuration 1 Swap %: 4.76 Reproduce %: 47.62 Select %: 47.62 ---------------- Configuration 2 Swap %: 33.33 Reproduce %: 33.33 Select %: 33.33 ---------------- Configuration 3 Swap %: 61.26 Reproduce %: 19.37 Select %: 19.37 ---------------- Configuration 4 Swap %: 83.33 Reproduce %: 8.33 Select %: 8.33 ----------------
Running simulations for 300 steps.
model11 = run_model(model_class=RSPMobility, params=model_params11, steps=300)
model12 = run_model(model_class=RSPMobility, params=model_params12, steps=300)
model13 = run_model(model_class=RSPMobility, params=model_params13, steps=300)
model14 = run_model(model_class=RSPMobility, params=model_params14, steps=300)
Visualizing the grid after 300 steps:
fig, axs = plt.subplots(2, 2, figsize=(25,15))
draw_grid(
model11.grid,
agent_portrayal,
title=build_title(params=model_params11, steps=300),
ax=axs[0][0]
)
draw_grid(
model12.grid,
agent_portrayal,
title=build_title(params=model_params12, steps=300),
ax=axs[0][1]
)
draw_grid(
model13.grid,
agent_portrayal,
title=build_title(params=model_params13, steps=300),
ax=axs[1][0]
)
draw_grid(
model14.grid,
agent_portrayal,
title=build_title(params=model_params14, steps=300),
ax=axs[1][1]
)
for ax in axs.flat:
ax.set_aspect('equal')
As in the model with simultaneous activation of agents, species self-arrange by forming patterns of moving spirals. Increasing the swap exponent (i.e. the mobility of agents), the spiral structures grow in size.
Visualizing the evolution of the system:
fig, axs = plt.subplots(2, 2, figsize=(25,10))
line_plot(
model11,
title=build_title(params=model_params11),
ax=axs[0][0]
)
line_plot(
model12,
title=build_title(params=model_params12),
ax=axs[0][1]
)
line_plot(
model13,
title=build_title(params=model_params13),
ax=axs[1][0]
)
line_plot(
model14,
title=build_title(params=model_params14),
ax=axs[1][1]
)
Let's define the parameters to perform a simulation of a system with 3 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the three species and empty sites in equal proportions (init0=0.25, init1=0.25, init2=0.25, init3=0.25), and setting the swap exponent to 1 and the reproduce and select exponents to -0.7.
model_params15 = {
"height": 100,
"width": 100,
"n_species": 3,
"init0": 0.25,
"init1": 0.25,
"init2": 0.25,
"init3": 0.25,
"init4": 0,
"init5": 0,
"swap_exp": 1,
"reproduce_exp": -0.7,
"select_exp": -0.7,
"color_map": color_map_mobility
}
This corresponds to the following percentages of events at each step:
compute_mobility_perc(model_params15)
Swap %: 96.16 Reproduce %: 1.92 Select %: 1.92
Running a simulation for 2000 steps.
model15 = run_model(model_class=RSPMobility, params=model_params15, steps=2000)
Visualizing the evolution of the system:
line_plot(
model15,
title=build_title(params=model_params15)
)
When mobility is high, biodiversity is jeopardized and lost: the system adopts a uniform state where only one species is present, while the others have died out. Which species remains is subject to a random process, all species having equal chances to survive in the model.
Let's define the parameters to perform a simulation of a system with 4 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the five species and empty sites in equal proportions and setting swap, reproduce and select exponents to 0.
model_params16 = {
"height": 150,
"width": 150,
"n_species": 4,
"init0": 0.2,
"init1": 0.2,
"init2": 0.2,
"init3": 0.2,
"init4": 0.2,
"init5": 0,
"swap_exp": 0,
"reproduce_exp": 0,
"select_exp": 0,
"color_map": {
0: 'black',
1: 'red',
2: 'purple',
3: 'gold',
4: 'green'
}
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPMobility,
model_params16,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the evolution of the species proportions during 300 steps.
model16 = run_model(model_class=RSPMobility, params=model_params16, steps=300)
draw_grid(
model16.grid,
agent_portrayal,
title=build_title(params=model_params16, steps=300)
)
line_plot(
model16,
title=build_title(params=model_params16)
)
With 4 species, the system dynamics stabilize with non competing species forming islands.
Let's define the parameters to perform a simulation of a system with 5 species, on a square lattice with $N=150\cdot 150$ sites, initialized with the sites being randomly assigned to each of the five species and empty sites in equal proportions and setting swap, reproduce and select exponents to 0.
model_params17 = {
"height": 150,
"width": 150,
"n_species": 5,
"init0": 1/6,
"init1": 1/6,
"init2": 1/6,
"init3": 1/6,
"init4": 1/6,
"init5": 1/6,
"swap_exp": 0,
"reproduce_exp": 0,
"select_exp": 0,
"color_map": {
0: 'black',
1: 'red',
2: 'purple',
3: 'gold',
4: 'green',
5: 'blue',
}
}
Creating the widget to interactively simulate the system:
page = JupyterViz(
RSPMobility,
model_params17,
measures=[line_plot],
agent_portrayal=agent_portrayal,
)
page
Running a simulation and plotting the evolution of the species proportions during 300 steps.
model17 = run_model(model_class=RSPMobility, params=model_params17, steps=300)
draw_grid(
model17.grid,
agent_portrayal,
title=build_title(params=model_params17, steps=300)
)
line_plot(
model17,
title=build_title(params=model_params17)
)
Also with 5 species, entangled rotating spiral waves emerge.
4. References¶
- [1] Frean, Marcus, and Edward R. Abraham. "Rock–scissors–paper and the survival of the weakest." Proceedings of the Royal Society of London. Series B: Biological Sciences 268.1474 (2001): 1323-1327.
- [2] Reichenbach, Tobias, Mauro Mobilia, and Erwin Frey. "Mobility promotes and jeopardizes biodiversity in rock–paper–scissors games." Nature 448.7157 (2007): 1046-1049.